Dipolar cross-relaxation in a two-spin system

Ilya Kuprov,
University of Oxford.
kuprov@physchem.ox.ac.uk

This is a simple program for calculation of diamagnetic dipolar relaxation and cross-relaxation rates between user-specified nuclear spin orders in a two-spin system. The internuclear distance is assumed to be fixed (rigid molecule). Specify the source and the target spin order in the program body and evaluate the entire worksheet. Pencil evaluation of these quantities is a rather long story, the best places to look are:

M. Goldman, Formal theory of spin-lattice relaxation, J. Magn. Res., 149, 160-187 (2001).
A.G. Palmer III, in: Protein NMR spectroscopy: principles and practice, Academic Press, 1996.
M. Mehring, V.A. Weberruss, Object-Oriented Magnetic Resonance, Academic Press, 2001.

The calculation uses group theoretical properties of the dipolar interaction Hamiltonian, expanding it into a sum of irreducible tensor operators, which are rotated using Wigner matrices. The Wigner matrix elements are then autocorrelated in the Bloch-Redfield-Wangsness theory treatment.

General setup:

Needs["LinearAlgebra`MatrixManipulation`"]; 
      Off[General::spell1, ClebschGordan::phy];
Id=(({{1, 0}, {0, 1}})); (Subscript[σ, x])=(({{0, 1/2}, {1/2, 0}})); (Subscript[σ, y])=(({{0, -I/2}, {I/2, 0}})); (Subscript[σ, z])=(({{1/2, 0}, {0, -1/2}}));(Subscript[T, 1, -1])=(({{0, 0}, {1/(Sqrt[2]), 0}}));(Subscript[T, 1, 0])=(({{1/2, 0}, {0, -1/2}}));(Subscript[T, 1, 1])=-(({{0, 1/(Sqrt[2])}, {0, 0}}));
Kron[a_, b_] := BlockMatrix[Outer[Times, a, b]];Comm[A_, B_] := A . B - B . A;Scal[A_, B_] := Tr[A . Conjugate[Transpose[B]]];
T[k_, 
          q_] := (((-1))^q)  (Sqrt[2 k + 1]) (∑_(q1 = -1)^1)(∑_(q2 = -1)^1)(ThreeJSymbol[{1, q1}, {1, 
                      q2}, {k, -q}] Kron[(Subscript[T, 1, q1]), 
                    (Subscript[T, 1, q2])]);
(Subscript[H, x])=Kron[(Subscript[σ, x]), Id];(Subscript[H, y]) = Kron[(Subscript[σ, y]), Id];(Subscript[H, z]) = Kron[(Subscript[σ, z]), Id];(Subscript[H, plus]) = (Subscript[H, x]) + I (Subscript[H, y]);(Subscript[H, minus]) = (Subscript[H, x]) - I (Subscript[H, y]);
(Subscript[F, x]) = Kron[Id, (Subscript[σ, x])];(Subscript[F, y]) = Kron[Id, (Subscript[σ, y])];(Subscript[F, z]) = Kron[Id, (Subscript[σ, z])];(Subscript[F, plus])=(Subscript[F, x]) + I (Subscript[F, y]);(Subscript[F, minus]) = (Subscript[F, x]) - I (Subscript[F, y]);
(Subscript[DQ, x]) = 
      (1/(Sqrt[2])) (2 (Subscript[H, x]) . (Subscript[F, x]) - 2 (Subscript[H, y]) . (Subscript[F, y])); 
    (Subscript[DQ, y]) = (1/(Sqrt[2])) (2 (Subscript[H, y]) . (Subscript[F, x]) + 2 (Subscript[H, x]) . (Subscript[F, y])); 
    (Subscript[ZQ, x]) = (1/(Sqrt[2])) (2 (Subscript[H, x]) . (Subscript[F, x]) + 2 (Subscript[H, y]) . (Subscript[F, y])); 
    (Subscript[ZQ, y]) = (1/(Sqrt[2])) (2 (Subscript[H, y]) . (Subscript[F, x]) - 
            2 (Subscript[H, x]) . 
                (Subscript[F, y]));
$Assumptions = τc > 
            0 && τ ∈ Reals && 
          t > 0 && (Subscript[ω, 1]) ∈ Reals && (Subscript[ω, 2]) ∈ 
            Reals && (Subscript[A, dip]) > 0 && ℏ > 0 && 
          Jc ∈ Reals;

Shorthands:

\(\(\(\ \)\(ω\_1 = \(-γ\_1\)\ B\_0, \ \ \ \ \ \ \ \ \ \ \ \ ω\_2 = \(-γ\_2\)\ B\_0, \ \ \ \ \ \ \ \ \ \ \ \ A\_dip = \(γ\_1\ γ\_2\ ℏ\ μ\_0\)\/\(4\ π\ r\^3\)\)\)\)

\(Source\ and\ target\ spin\ orders\ in\ the\ laboratory\ \(\(frame\)\(:\)\)\)

Source = (Subscript[H, z]);
Target = (Subscript[F, z]);

Bloch - Redfield - 
          Wangsness theory processing:

Hst = (Subscript[ω, 1]) (Subscript[H, z]) + (Subscript[ω, 2]) (Subscript[F, z]) // 
          Simplify;
Hdn[t_]=Total[((({{(Subscript[D, -2, -2])[
                              t], (Subscript[D, -2, -1])[
                              t], (Subscript[D, -2, 0])[
                              t], (Subscript[D, -2, 1])[
                              t], (Subscript[D, -2, 2])[t]}, {(Subscript[D, -1, -2])[
                              t], (Subscript[D, -1, -1])[
                              t], (Subscript[D, -1, 0])[
                              t], (Subscript[D, -1, 1])[
                              t], (Subscript[D, -1, 2])[t]}, {(Subscript[D, 0, -2])[t], (Subscript[D, 0, -1])[
                              t], (Subscript[D, 0, 0])[t], (Subscript[D, 0, 1])[
                              t], (Subscript[D, 0, 2])[t]}, {(Subscript[D, 1, -2])[t], (Subscript[D, 1, -1])[
                              t], (Subscript[D, 1, 0])[t], (Subscript[D, 1, 1])[
                              t], (Subscript[D, 1, 2])[t]}, {(Subscript[D, 2, -2])[t], (Subscript[D, 2, -1])[
                              t], (Subscript[D, 2, 0])[t], (Subscript[D, 2, 1])[
                              t], (Subscript[D, 2, 2])[t]}})).(({{0}, {0}, {(Sqrt[6])}, {0}, {0}})))*(({{T[2, 2]}, {T[2, 1]}, {T[2, 0]}, {T[2, -1]}, {T[2, -2]}})),2]//Simplify;
Dcomm = -Comm[Hdn[t], 
              Comm[Conjugate[
                  Transpose[
                    MatrixExp[I Hst τ] . Hdn[t + τ] . 
                      MatrixExp[-I Hst τ]]], 
                Source]] // 
          Simplify;
Rate = Scal[Dcomm, Target]/(Sqrt[Scal[Target, Target] Scal[Source, Source]]) // Simplify // 
          ExpandAll;
For[a = -2, 
      a ≤ 2, 
      For[b = -2, b ≤ 2, 
        For[c = -2, c ≤ 2, 
          For[d = -2, 
            d ≤ 2, 
Rate = 
              ReplaceAll[
                Rate, {(Subscript[D, a, b])[t] Conjugate[
                        (Subscript[D, c, d])[t + τ]] -> 
                    KroneckerDelta[a, c] KroneckerDelta[b, 
                        d] (1/5) (E^(-τ/τc))}]; 
d++;]; c++;]; b++;]; a++;]; 
    Clear[a, b, c, d];
 Rate = 
          (Subscript[A, dip]^2) Integrate[Rate, {τ, 0, ∞}, 
              GenerateConditions -> 
                False] ;Rate = ComplexExpand[Rate] // Re // 
            FullSimplify ;

Source to target cross - relaxation rate:

Rate